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ABSTRACT 

The decomposition of an image into a linear combination of digitised basis functions is an ev- 
eryday task in astronomy. A general method is presented for performing such a decomposition 
optimally into an arbitrary set of digitised basis functions, which may be linearly dependent, 
non-orthogonal and incomplete. It is shown that such circumstances may result even from the 
digitisation of continuous basis functions that are orthogonal and complete. In particular, digi- 
tised shapelet basis functions are investigated and are shown to suffer from such difficulties. 
As a result the standard method of performing shapelet analysis produces unnecessarily inac- 
curate decompositions. The optimal method presented here is shown to yield more accurate 
decompositions in all cases. 

Key words: methods: data analysis - techniques: image processing. 
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q . 1 INTRODUCTION 

h ; 

| The linear decomposition of an astronomical dataset into a series of 
C^J . basis functions is a fundamental task in many areas of astrophysics 
!L* ' and cosmology. Indeed, in the analysis of one-dimensional spec- 
. ,_( I tra, two-dimensional images or higher-dimensional datasets, one 
■ is often faced with the problem of determining the coefficients in 
•_h ' a linear expansion of the dataset in some set of (sampled) basis 
functions. To illustrate our discussion, we will focus here on the 
important specific example of decomposing a two-dimensional as- 
tronomical image, although the general approach that we advocate 
will be applicable to datasets of arbitrary dimensionality. 

The description of a digitised image as a linear combination 
of a set of sampled basis functions is an everyday problem for as- 
tronomers. For example, one often describes an image in terms of 
a set of orthogonal Fourier modes by calculating the coefficients in 
the expansion using a discrete Fourier transform. Significant atten- 
tion has also been given to representing an image as a linear com- 
bination of discrete wavelets basis functions, both orthogonal and 
non-orthogonal (see, for instance, Hobson et al. 1999; Sanz et al 
1999a;b; Tenorio et al. 1999). More recently, several authors have 
investigated the use of Gaussian-Hermite modes (or shapelets) in 
representing images of galaxies (Refregier 2003; Refregier & Ba- 
con 2003). Although such image decompostion is commonplace, 
there exist a number of subtleties in the procedure that are not 
widely appreciated within the astronomical community. These in- 
clude the effect of digitisation, or sampling, on familiar notions 
from the analytic theory of continuous functions, in particular the 
concepts of completeness and linear independence of sampled basis 
functions, and the use of non-orthogonal bases. 

In this paper, we therefore discuss a general mathematical 
framework for the linear decomposition of an astronomical im- 



age into a set of arbitrary sampled functions (or modes). These 
may, in general, form a 'basis' that is non-orthogonal and either 
under-complete, perfectly-complete or over-complete. In the under- 
complete case, the modes do not support all the degrees of free- 
dom in an arbitrary image. By perfectly-complete we mean that the 
modes supports all the degrees of freedom without linear depen- 
dence between them, whereas for an over-complete mode set there 
are additional modes relative to the complete case and hence lin- 
ear dependence in the set. In particular, we show how to quantify 
those degrees of freedom in an image that can be supported by any 
particular mode set, and we explicitly consider how to obtain the 
optimal set of mode coefficients for representing the image. These 
methods present an opportunity for the optimisation of sampling 
within modal analysis, ensuring that computational operations are 
reduced to a minimum. 

The structure of the paper is as follows. In Section|2| we out- 
line the modal decomposition problem and present a general pre- 
scription for obtaining optimal results, which is based on the sin- 
gular value decomposition of the matrix of digitised mode vectors. 
In Section [5] we illustrate our general approach by investigating 
the decomposition of simple one-dimensional function into Gauss- 
Hermite (shapelet) modes. This investigation is pursued further in 
Section|4| where we consider the decomposition of images of Hub- 
ble Deep Field galaxies into two-dimensional Gauss-Laguerre (po- 
lar shapelet) modes. Finally, our conclusions are presented in sec- 
tion|5] 



2 THE MODAL DECOMPOSITION PROBLEM 

Consider a 2-D digitised image consisting of N pixels, which we 
denote by the ^-dimensional column vector d. Our goal is to rep- 
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resent the image, as accurately as possible, in terms of a set of 
M (two-dimensional) modes {e^} (k = 1,2,... ,M) defined at the 
same points at which d is sampled; thus each is also a column 
vector of length N. It is convenient to combine the mode vectors 
into the single N xM mode matrix defined by 

E= [e\ ■■■e k ---e M \- 

Our objective is thus to determine a suitable coefficient vector a 
that minimises the residual 

e=\Ea-d\. (1) 

Note that we make no requirement for the mode vectors to be mu- 
tually orthogonal, linearly independent, or of unit length. 

Several cases may arise in this problem. There may exist a 
unique coefficient vector a that minimises Q; this corresponds to 
the quantity e 2 possessing a single well-defined multi-dimensional 
quadratic minimum in the space of coefficients a. Surfaces of 
constant £ 2 are given by the Hermitian form (a — a)' R(a — a) = 
constant, where R =E'E (which is called the Gram matrix). The 
eigenvectors of R determine the principal directions of this multi- 
dimensional ellipsoid, and the extent along each principal direction 
is inversely proportional to the corresponding eigenvalue. If there 
exist approximate degeneracies in the coefficient-space, some of 
the eigenvalues become very small and so the multi-dimensional 
ellipsoid is considerably elongated in these directions; this leads to 
a wide range of coefficients vectors a for which e takes its mini- 
mum value to within the numerical precision. In the limiting case 
of exact degeneracies, some eigenvalues are identically zero lead- 
ing to an infinite extent for the 'ellipsoid' along the corresponding 
principal directions. In the subspace spanned by these directions the 
value of e takes its minimum value precisely. In either of the last 
two cases, to the numerical precision, there exist an infinite num- 
ber of possible coefficient vectors a that minimise 0. The value 
of £ m ; n is also of central importance. Clearly, if e m j n = 0, then the 
modal expansion Ed provides an exact representation of the origi- 
nal image, whereas if > it represents only an approximation 
tod. 

Which of the above cases occurs is dependent on the com- 
pleteness and linear dependence of the mode set ej., and on the 
particular image d being decomposed. Clearly, if the M < N then 
the modes are 'under-complete', and so cannot represent an arbi- 
trary image, but only an approximation to it. This may also occur 
when M >N, however, since the modes may still not form a com- 
plete basis as a result of linear dependence between them. Nev- 
ertheless, even with an under-complete set, the particular image d 
under analysis may lie in the subspace spanned by the modes, and 
so may be represented exactly. When M >N cases arise where the 
modes form a complete basis, in which any image can be repre- 
sented exactly. In this case, if M — N the modes must linearly in- 
dependent and form a 'perfectly-complete' basis. If M > N, some 
linear dependence must exist between the modes and they form an 
'over-complete' basis. 

In general, the degree of completeness of the chosen mode 
set may not be known at the outset. Fortunately, a straightforward 
technique exists for determining the degree of completeness of the 
mode set, and simultaneously determining an 'appropriate' coef- 
ficient vector a that minimises Q. This is achieved by perform- 
ing the singular- value decomposition (S VD) of the mode matrix E, 
which we now discuss. 



2.1 Singular value decomposition 

The SVD of the N x M mode matrix E (which may, in general, 
contain complex-valued entries) may be written (see, for example, 
Golub & Van Loan 1992) 

E = UI.V\ (2) 

where U is unitary matrix of dimensions N x N, V is a unitary ma- 
trix of dimensions M x M, and £ is a A' x M matrix (the same di- 
mensions as E) that is diagonal in the sense that = o, for i < p, 
where p = min[M,N], and zero otherwise. The coefficients a, are 
the singular values of the matrix E. 

As is well-known, the number r (say) of non-zero singular val- 
ues is equal to the rank of E, which in turn is the dimensionality of 
the image subspace that is spanned by the M modes {e^}; clearly it 
must be the case that r <M and r < N. The dimensionality of the 
nullspace of E (the subspace of vectors n in the coefficient space for 
which En = 0) is given by M — r. It is thus a simple matter to deter- 
mine the completeness of the mode set contained in E as follows: 
(i) if r < N the modes are under-complete; (ii) if r = N and M = N 
the modes are perfectly complete; and (iii) if r — N and M > N the 
modes are over-complete. Moreover, it is straightforward to show 
that the columns of U corresponding to non-zero singular values 
constitute an orthonormal basis for the range of E, and the columns 
of V that do not correspond to a non-zero singular value form an or- 
thonormal basis for the nullspace. In practical numerical problems, 
it is often the case that none of the singular values a, are iden- 
tically zero. Instead, one usually sets to zero those singular values 
for which |o", |/|o"i | < r|, where r) is some small factor (for example, 
10~ 5 in single precision arithmetic) and Oi is the first (and largest) 
singular value. 

Once the SVD J2j has been calculated, it is also straightfor- 
ward to obtain an 'appropriate' coefficient vector a that minimises 
@. As we outline below, in all cases one should calculate the coef- 
ficient vector using 

a = VZ f UU, (3) 

where the M xN matrix denoted by £ is constructed by taking the 
transpose of £ in J2j and replacing each non-zero singular value a,- 
by I/O;. 

It is clear that, with the above construction, ££"* is an N x N 
diagonal matrix with diagonal entries that equal unity for those val- 
ues of j for which Oj ^ 0, and zero otherwise. Using this result it is 
straightforward to show that ^3} does indeed minimise the residual 
Q. Modifying slightly the argument of Press et al. (1994), suppose 
we were to add to a some arbitrary vector a' = a[ +a' 2 , where a' 2 is 
the part of the vector that lies in the nullspace of E (if one exists) 
and a'j is the part that lies in the complement to the nullspace. This 
would result in the addition of the vector d' = Ea[ to Ea d. We 
would then have 

\Ea-d+d'\ = | (l/IE t I/ 1 " - I)d + d! 

= [£/[(££ + -/)£/ t rf+£/ t ]rf / | 

= \{ZZ-I)uU + U ] 'd'l (4) 

where in the last line we have made use of the fact that the length 
of a vector is left unchanged under the action of the unitary matrix 
U. Now, the jth component of the vector (T.L — I)U^d will only 
be non-zero when Cj = 0. However, the jth element of the vector 
U^d' is non-zero only if 0^ 7^ 0, since d' lies in the range of E. 
Thus, as these two terms only contribute to ^4) for two disjoint sets 
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of y'-values, its minimum value, as a 1 is varied, occurs when d' = 0; 
this requires a\ = 0. 

If the image d lies in the subspace spanned by the M modes 
{e^}, then the minimum value of the residual Q is zero, and the 
image is represented exactly by Ea. Clearly, in order to represent an 
arbitrary image we require r = N (and so M > N). If d does not lie in 
the range of E, then Ea yields only an approximate representation 
of the image (for this to occur, the modes must be under-complete; 
this will always be true if M < N, but may also occur when M > N). 

In either case, we see from the argument above that, if E does 
not possess a null space, the coefficient vector J5) is unique in 
minimising the residual Q. The condition for this to occur is that 
r = M, indicating that the M modes are linearly independent (and 
hence requiring M < N).\f E does possess a nullspace, however, 
any vector a' 2 in this nullspace can be added to J3j, without chang- 
ing the value of the residual. Thus there are an infinite number of 
coefficient vectors that minimise Q. The condition for this to oc- 
cur is that r < M, indicating linear dependence between the modes 
(which may occur for M <N and M > N). 

In the case where E possesses a nullspace, from the infinite 
number of coefficient vectors that minimise the residual Q, the 
vector a in J3J is that which contains no contribution from the 
nullspace. An equivalent statement is that J3J is the vector of short- 
est length that minimises Q. Consider again adding some arbitrary 
vector a' to J3j. The length of the resulting vector is 

\a+a'\ = \VL f uU+a'\ 

= |F(zVrf + vV)| 

where in the last line we again make use of the fact that the length 
of a vector is left unchanged under the action of a unitary matrix. 
The jth component of the vector LU'd will only be non-zero when 
Cj 7^ 0, whereas the jth element of the vector V a! is non-zero only 
if Oj = 0. Thus the minimum length vector has a' = 0. 

2.2 Dual modes 

Although the appropriate coefficient vector may always be obtained 
straightforwardly from it is conceptually appealing to consider 
coefficient as the scalar product of the image vector d with some 
new mode vector that is dual to the original mode e^. In other 
words, it is often convenient to rewrite as 

a=E f d, (5) 

where E is the N xM dual mode matrix, which contains the M dual 
mode vectors as its columns, i.e. 

E= [S\ ■■■e k ---e M \. 

From J3}, we see immediately that the dual mode matrix is given in 
terms of the SVD of the original mode matrix by 

E = Vtv\ 

Thus, once a mode set has been defined, the dual mode set can be 
calculated directly by performing a SVD, without reference to any 
image. The coefficients for any particular image are then quickly 
obtained using 

It should be noted that, in general, 

EE 1 ^I N ^EE\ 

where we now start to place subscripts on identity matrices to em- 
phasise their dimensionality. It is only in the case where the original 



mode set is complete that EE = Iff = EE ' , from which it is a sim- 
ple matter to verify that EE = /# — EE^ . This, of course, corre- 
sponds to the case in which an arbitrary image d can be represented 
exactly, since then 

Ea=EE ! d = d. 

Otherwise only an approximation to the original image is possible. 
It is also worth pointing out that, in general, 

E*E ^I M ¥= E f E. 

In other words, the dual modes do not necessarily obey the stan- 
dard orthogonality condition etei = 8*/ = e' k ei with the original 
mode set. This orthogonality condition is only satisfied if the orig- 
inal mode set is linearly-independent (for which r = M and hence 
M < N). In this case, E^E =Im = E^E, from which one quickly 
finds that E^E = Im —E E, which recovers the standard result that 
the dual modes form the reciprocal basis of the original mode set. 
This corresponds to the case where {5J is unique in minimising the 
residual Q. If, in addition, we have M — N and so the basis per- 
fectly complete, E is invertible and it follows that E = E~ l . We 
reiterate, however, that all possible cases are automatically accom- 
modated using SVD. 

2.3 Orthogonal mode sets 

So far, we have imposed no restrictions on the orthogonality, nor- 
malisation or the number of members of our original mode set {e^}. 
It is worth investigating, however, the simplifications that occur in 
the case where the modes are mutually orthogonal and each mode 
is unique (so M < N). Hence, we have E'E = Im- 

In this case, the mode set is automatically linearly- 
independent and so the coefficient vector {3} or |5) is unique in min- 
imising the residual. Moreover, the dual modes and original modes 
obey the standard orthogonality condition E^E = Im =E*E.\& the 
case of orthogonal modes, it is clear that, in addition, this condition 
is satisfied if the dual modes are given by 



Thus, for orthogonal modes (even with M < N), we see that each 
dual is simply proportional to the original mode. In the event that 
the original modes are normalised to unit length, this proportion- 
ality becomes an equality, and the dual modes are identical to the 
original modes, so that E = E and hence 

a = EU. (6) 

This last result is, of course, the reason for the common practice, 
in the case of orfhonormal modes, of calculating the kth mode co- 
efficient by simply evaluating the scalar product of the image d 
with the kth mode e^. Of course, when M = N, a set of orthogonal 
modes also forms a perfectly complete basis. 

2.4 Practical considerations 

Throughout our discussion, we have adopted the practical approach 
of working with finite dimensional vectors, rather than continous 
functions. Thus, both the image d and the modes are considered 
simply as W-dimensional vectors. This allows the direct application 
of our approach to a wide variety problems in the decomposition of 
digital astronomical images. The image pixel values d, correspond 
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to samples of the underlying continuous distribution d(x) at a par- 
ticular set of points {*, }. In practice, it is most often the case that 
sample points are regularly spaced, but this restriction is not neces- 
sary. The entire approach can just as easily be applied to the case 
where the pixels values correspond to samples that are not regularly 
spaced. In an analogous manner, the elements of each mode vector 
are the sample values from some underlying continuous mode 
function fkix) at the same set of points {.*,■}. 

It is clear that the positions of the sample points will have a 
profound effect on whether the linear-independence, orthogonality, 
or other properties of the original set of continuous mode func- 
tions {fk(x)} are inherited by the mode vectors {e^}- Consider, for 
example, a set of continuous mode functions {fk(x)} that are or- 
thonormal over some continuous domain D , so that 



fk*{x)fl{x) d* = 5* 



Even in the case where the sample points {x,} are evenly spaced, 
the spacing between them and the extent of the domain £> that is 
sampled are important in determining whether the corresponding 
mode vectors {e^} are orthonormal. Indeed, as we shall see below, 
in some applications it is common for the corresponding mode vec- 
tors not to be orthonormal. If this is the case, then it is is no longer 
true that the dual mode vectors are identical to the original mode 
vectors, and so it is incorrect to calculate the coefficient vector a us- 
ing J6j, i.e. by taking the scalar product of the datarf with the mode 
vectors. This will lead to a modal decomposition of the image that 
is unnecessarily inaccurate. Instead, one must return to using the 
general result J5J, where the mode coefficients are calculated by 
taking the scalar product of the data with the corresponding dual 
mode vectors. 

Finally, a crucial point is that the converse may also occur. For 
example, if a set of continous mode functions are not orthogonal, 
by sampling them at a particular set of (non-uniform) points one 
can arrive at a set of mode vectors that are orthogonal, in which 
case the duals are trivially obtained. 



3 ONE-DIMENSIONAL ILLUSTRATIONS 

Although the main focus of this paper is the modal decompostion 
of two-dimensional astronomical images, it will be informative first 
to illustrate the discussion given above by investigating the decom- 
postion of one-dimensional sampled functions (which may be con- 
sidered as cuts through some image). In particular, we will focus 
our numerical investigations in this section on decompostions using 
one-dimensional Gaussian-Hermite (or shapelet) mode functions, 
as recently advocated by Refregier (2003) for representing images 
of galaxies. 

In this case, the continuous mode functions are given by 

A(r,p) = p- 1/2 te(P -1 *), 

for k = 0, 1,2, . . . The quantity P is some fixed 'characteristic 
scale' for the set of mode functions and the function §k( u ) reads 

<Sfk{u) = {2 k n^ 2 k\)-^ 2 H k {u)exp{-u 2 /2), 

where Hk(u) is a Hermite polynomial of order k. From these ex- 
pressions it is clear that the mode functions are real and that P is 
equal to the dispersion 'o"' of the k = Gaussian mode function. 
It is also straightforward to show that the mode functions are or- 
thonormal over the domain — o° < x < °°, i.e. 
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Figure 1. The Gaussian-Hermite continuous mode functions ft (.v;|3) for 
it = 0,1,..., 5 and |3 = 12. 



We note that, for convenience, the mode functions are centred at 
x — 0. It is, of course, possible to centre them about some arbitrary 
point x = x c , but in most cases the coordinate system is chosen so 
that the centroid of the galaxy image (say) is located at the origin. 
The first six mode functions are shown in Fig.Qfor p = 12. 

Let us now consider the decomposition of a one-dimensional 
pixelised 'image' into shapelets. In any such decomposition, one is 
free to choose both the characteristic scale P of the modes set and 
the total number M of modes used (note that M = k m?LX + 1). As 
shown by Refregier (2003), such a set of shapelet modes is suitable 
for describing features in an image on scales between the two limits 
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Thus, if the image has features on scales ranging from 9 m j n (e.g. the 
pixel size or the size of the point spread function) to max (e.g. the 
size of the image or the overall extent of the structure it contains), 
then a good choice of p and M is 
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/*(.r;p)/K*;P)<k = S fc 



(7) 



For each of the pixelised test images we consider below, we 
take the coordinate positions of the N (regularly-spaced) sample 
points to be 

{*,-} = {4 (iV-1),..., -1,0,1,... ,|(iV-l)}, 

thereby giving an image of spatial extent of N, with a pixel size of 
unity. In each case, we take iV = 81. Taking into account that our 
test images contain structure on the order of the image size, but do 
not contain structure on scales less than a few pixels, acceptable 
choices for the scale parameter and the number of modes to use are 
P ~ 10 and M ~ 20. For convenience, we fix M = 20, and investi- 
gate mode sets for which P = 1 - 20 in steps of unity. In each case, 
the elements of the corresponding pixelised mode vectors are easily 
calculated from (ej); = P). 

In Fig.|2|we plot the elements of the first six mode vectors for 
p = 12. We see that, even for this modest value of P, the support of 
the mode vectors for k > 2 exceeds the length of the image. Indeed, 
as k increases the support of the mode vectors increases according 
to {5}, and so this effect becomes extremely pronounced for the 
high-A; modes. As a result, we would expect the mode set does not 
satisfy the discretized form of the orthogonality condition Q. We 
will show below that this is indeed the case. 

It is also of interest to investigate the structure of the mode 
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Figure 2. The Gaussian-Hermite mode vectors for k = 0, 1,... ,5 and 
p- 12. 
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Figure 3. The i max = 19 Gaussian-Hermite mode vector e\<) for p 
0,1,2,3,6,10,15. 



vectors for different values of (3. As an illustration, in Fig.[3]we plot 
the highest-A: mode vector in the set (fc max = 19) for six different 
values of the scale parameter p. In particular, we note that the mode 
vector is severely undersampled for p = 1 and P = 2, in which case 
the mode set may again fail to be orthogonal. Moreover, since the 
highest-A: mode in the set has the largest spatial extent, we see that, 
in fact, for P > 6 the support of the mode set exceeds the length of 
the image, once again suggesting a loss of orthogonality. 
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Figure 4. The rank r of the mode matrix E as a function of the scale param- 
eter P for M = 20 Gaussian-Hermite mode vectors. 



3.1 Properties of the mode vectors 

Since we have chosen M <N, the mode vectors clearly cannot form 
a complete basis for the image space. Nevertheless, it is of interest 
to investigate formally the linear-dependence and orthogonality of 
the mode vectors e$ [k = 0, 1, ... ,M— 1) for various values of the 
scale parameter p. To this end, for each value of p under considera- 
tion, we construct the corresponding mode matrix E and then deter- 
mine its rank r by calculating its SVD and counting the number of 
non-zero singular values. If the mode vectors are linearly indepen- 
dent, one would expect r = M. In Fig.|4] we plot the value of r as 
a function of the parameter p. Since the number of mode vectors is 
M = 20, we see that the mode vectors form a linearly-independent 
(r = 20) set only if the scale parameter lies in the range P = 2—11. 
For smaller values of P, the undersampling of the underlying con- 
tinuous mode functions leads to linear dependence in the resulting 
mode set. For large P, linear dependence occurs since the support 
of (some of) the mode vectors exceeds the length of the image. 

To investigate whether the mode vectors are orthonormal we 
calculate the Gram matrix R=E^E for each value of P under con- 
sideration. Clearly, for an orthonormal set of mode vectors would 
expect R = Im- In Fig [5] (top), for each value of P, we plot the 
logarithms of the largest and smallest diagonal elements of R . For 
mode vectors of unit length, we would expect all the diagonal ele- 
ments to equal unity, and so the two plotted curves should coincide 
at log (/?,-,-) = 0. We see that this occurs only if the scale parameter 
lies in the range P = 3-6. In the bottom panel of Fig. [5] we plot the 
value of the logarithm of the largest off-diagonal element of R as a 
function of p. For an orthogonal set of mode vectors, one would ex- 
pect all the off-diagonal elements of R to be zero in the ideal case, 
or below the machine precision in a numerical implementation. We 
see that this is only the case if P = 3-4, although the off-diagonal 
elements remain reasonably small for P = 5-6. For other values of 
p, however, we see that there exist off diagonal elements with ab- 
solute values greater than ~ 0.1, which is a clear indication that 
the corresponding mode sets are not orthogonal. We also note that 
there exist values of p for which the modes are non-orthogonal, but 
still have full rank. 



3.2 Dual mode vectors 

Since for mode sets with p < 3 or p > 6 the mode vectors are not 
orthogonal, the corresponding dual mode vectors are not simply 
some multiple of the original modes. It is therefore of interest to 
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Figure 5. Top: the logarithm of the value of the largest (solid line) and 
smallest (dashed line) diagonal elements of the Gram matrix R = E^E as 
a function of the scale parameter p. Bottom: the logarithm of the largest 
off-diagonal element of R as a function of p. 



investigate the form of the dual modes in these cases. As an illus- 
tration, in Fig.[6|we plot the first six dual mode vectors for the case 
(3=12. These should be compared with the corresponding origi- 
nal mode vectors plotted in Fig. Q We see that the form of each 
dual is very different from the corresponding original mode vector. 
In particular, we note that the form of the duals is determined by 
the entire original mode set. As a result, the large spatial extent of 
the large-fe original mode vectors has an effect on the form of the 
low-fc dual mode vectors, even though the support of correspond- 
ing low-fc original mode vectors does not exceed the length of the 
image. Thus, for example, the k = and k = 1 dual vectors are 
markedly different from the correspondinig original mode vectors, 
even though the latter lie completely within the image length. 

Although clear from the discussion in Section l2~2l it is worth 
pointing out once more that the dual modes are not the basis set in 
terms of which an image is described. In the approach presented, 
the image is still represented as a linear combination of the origi- 
nal Gauss-Hermite (or shapelet) mode vectors. It is simply that the 
value of the coefficient in the linear combination is given by 
the scalar product of the image vector d with the Mi dual mode 
e^, rather than with the original mode vector e^. Thus, the advanta- 
geous properties of the Gauss-Hermite modes (or any other mode 
set under consideration) can still be used in the analysis of the re- 
sulting modal decomposition. From Fig|6| however, it is clear that 
the dual vectors may not possess the same localisation as the basis 
vector to which they relate. Consequently, for the k = mode for 
example, a feature towards the edge of the field will affect the co- 
efficients of the mode even though the original mode vector falls to 
zero there. 
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Figure 6. The Gaussian-Hermite dual mode vectors et for k = 0, 1, . . . , 5 
and p = 12. 



3.3 Decomposition of simple functions 

Now that we have discussed the properties of the mode vectors and 
their duals for various values of the scale parameter p, it is of in- 
terest to investigate the effect of taking proper account of linear de- 
pendence and non-orthogonal. To this end, we perform the modal 
decomposition of some simple one-dimensional test functions us- 
ing both the method advocated in Section l2~T1 and the standard ap- 
proach in which the coefficients in the decompostion are obtained 
simply by taking scalar products of the image vector with the origi- 
nal mode vectors (see e.g. Refregier 2003). The test functions con- 
sidered, although simple, are relevant to the decomposition of as- 
tronomical images. 

Uniform function 

We begin by considering the modal decomposition of a uniform 
signal. The issue of accommodating (nearly) uniform background 
emission is clearly relevant in an astronomical context. The modal 
decomposition of this function was calculated, using both the stan- 
dard method and the SVD method, for M = 20 modes and for inte- 
ger values of the scale parameter in the range P = 1-20. In Fig.0 
we plot the resulting decompositions for the standard method (left 
column) and the approach advocated here (right column) for some 
representative values of p. These values have been chosen to corre- 
spond to mode sets for which the mode vectors are linearly depen- 
dent and non-orthogonal, as demonstrated in Fig. [5] In Fig. [8] we 
plot the residual £ = \Ea — d\ of the decompostions as a function of 
(3 for the standard method (solid line) and the SVD method (dashed 
line). 

We see from Fig. that, as one would expect, the decompo- 
sitions produced by the two methods differ, in some cases signif- 
icantly. It was confirmed, however, that for mode sets with P = 3 
and P = 4, for which the mode vectors are linearly-independent and 
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Figure 7. The uniform function (dashed line) and its modal decomposi- 
tion (solid line) into M = 20 Gaussian-Hermite mode vectors with scale pa- 
rameter (3. The decomposition coefficients are obtained using the standard 
method (left column) and the SVD method (right column). 
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Figure 8. Decomposition residual £ = \Ea — d\ as a function of (5 for the 
modal decomposition of the uniform function using the standard method 
(solid line) and the SVD method (dashed line). 



orthonormal, the modal decompositions produced by the two meth- 
ods coincide to machine precision. In the cases illustrated in Fig.0 
we see that the SVD method provides a decomposition that is con- 
sistently better than that of obtained with the standard approach. 
In particular, we note that for a wide range of P-values, the SVD 
method produces a decomposition that is visually indistinguishable 
from the input function, whereas the standard approach exhibits 
a pronounced ringing for all values of (3. A quantitative descrip- 
tion of the improved decomposition quality is given by Fig. 00 As 
expected, the SVD approach always produces the smallest decom- 
position residuals, indeed by many orders of magnitude for P>7. 
Since the decomposition of an image is a linear process, the ability 
to describe accurately a uniform function enables better decompo- 
sitions of discrete objects in a uniform background, which is clearly 



important in astronomical applications. We also note that, for the 
SVD method, the decomposition residual is a monotonically de- 
creasing function of the scale parameter p. For the standard method, 
however, there is a shallow minimum in the residual at P = 7. It is 
interesting that, from Fig. [5] for this value of p the mode set is nei- 
ther linearly-independent nor orthogonal. 

Top-hat function 

As our second illustration, we consider a top-hat function of width 
equal to 20 units. Once again, this simple function is relevant to 
astronomical image analysis of, for example, saturated images of 
compact objects. The modal decomposition of this function was 
again calculated, using both the standard method and the SVD 
method, for integer values of the scale parameter in the range 
P = 1-20. In Fig.|5|we plot the resulting decompositions for the 
standard method (left column) and the SVD method (right column) 
for some representative values of p. The decomposition residuals 
are plotted in Fig. llOl as a function P for the standard method (solid 
line) and the method (dashed line). 

We see from Fig. [5] that, once again, the two methods pro- 
duce different decompositions for the values of P illustrated, and 
that the SVD method produces superior results. For both methods, 
however, the presence of only M = 20 modes leads to considerable 
ringing in the decompostions for almost all values of p. Neverthe- 
less, it is interesting that for P = 1 the standard method is unable to 
produce a reasonable reconstruction, whereas an excellent decom- 
position is obtained over a limited extent of the function using the 
SVD method, the extent corresponding to the span of the mode set. 
The quantitative comparison of the two methods given in Fig. 1101 
once again shows that, by design, the SVD method always pro- 
duces the most accurate decomposition. In this case, however, the 
difference between the two methods is not as great as it was for 
the decomposition of the uniform function. We also note that the 
residuals for both methods in this case are minimised at P = 3. At 
this point, the minimum residual is identical for the two methods, 
since the original mode set is linear-independent and orthonormal 
to machine precision. 

$-model profile 

Our final one-dimensional illustration is that of a P-model profile. 
Clusters of galaxies are often modelled (King 1972) as spherically- 
symmetric, with an electron density profile of the form 

-36/2 

/ r \ " 

1 



n e (r) =n 



r c 



where r c is the cluster core radius and b is a constant (this is usually 
denoted by P, but this has already been used in this paper to denote 
the scale parameter of the mode set). In this case, one easily finds 
that the projected thermal Sunyaev-Zerdovich profile of the cluster 
is given by 

-X 



f(x)=f 1 + 



(10) 



where x the angular separation on the sky and X = (3b — 1 ) /2. As- 
suming the standard value b — 2/3 gives X — 1/2, and we use this 
value for our test function, together with a core length x c = 3. 

The modal decomposition of this function was calculated, us- 
ing both the standard method and the SVD method, for integer val- 
ues of the scale parameter in the range P = 1-20. In Fig. ll ll we plot 





the resulting decompositions for the standard method (left column) 
and the SVD method (right column) for some representative values 
of p. The decomposition residuals are plotted in Fig. ll2l as a func- 
tion P for the standard method (solid line) and the SVD method 
(dashed line). 

From Fig. l 111 we see again that the SVD methods produces su- 
perior decompositions for the values of P illustrated. For both meth- 
ods, some ringing is observed in the decompositions, but this is 
much less pronounced for the SVD method. In particular, for large 
values of P, the standard method cannot accurately reproduce the 
narrow peak of the P-model, whereas this is achieved by the SVD 
method. We also note that, for P = 1, the SVD method reproduces 
the central sharp peak very accurately. The quantitative comparison 
of the two methods given in Fig. II ll clearlv shows the improvement 
obtained using the SVD method, which produces a residual signif- 
icantly below that obtained using the standard method for p > 8. 
This is also the value of the scale parameter for which both meth- 



ods yield their lowest residual, and the corresponding set of original 
mode vectors are close to linearly-dependent and non-orthogonal. 



4 IMAGE DECOMPOSTION 

We now consider the modal decomposition of real two-dimensional 
astronomical images. These images were obtained by first using 
SExtractor (Bertins & Arnouts 1996) on the Hubble Deep Fields 
(HDFs; Williams et al. 1996, 1998). The convolution mask and de- 
tection parameters were adapted from those used by Massey et al. 
(2003). In particular, in order to allow the recovery of faint ob- 
jects, we adopted a low signal-to-noise detection threshold, DE- 
TECT.THRESH, of 1.3. Objects with CLASS.STAR > 97 per cent 
were discarded, since we wish to analyse only galaxies. The im- 
age was then segmented into small square 'postage stamp' regions 
around the remaining galaxies. The size of the images were set to 
51x51 pixels. In this section, we investigate the decomposition of 
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the HDF galaxy images into a set of two-dimensional Gaussian- 
Laguerre modes. 



4.1 Gauss-Laguerre mode functions 

The two-dimensional Gauss-Laguerre continuous mode functions 
are the simultaneous eigenstates of energy and angular momentum 
for the two-dimensional isotropic quantum harmonic oscillator. In- 
terestingly, these mode functions are also the solutions in polar co- 
ordinates to the paraxial wave equation in optics (see e.g. Gold- 
smith 1998). The relationship between the standard forms of these 
mode functions used in the above contexts is discussed in detail in 
the Appendix. 

In the field of optical and quasi-optical systems, it is custom- 
ary (Murphy et al. 1996; Withington et al. 2000) to define the set of 
mode functions as 

1/2 



<(r,e) = 



(2-8o m )p! 
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x H/2 e -. t /2 L M We i m6] (H) 



where x = r 2 /(3 2 and Lp (x) are the standard associated Laguerre 
polynomials. The functions \|/™ are conventionally called the Gaus- 
sian beam pm-modes or simply the pm-modes (see e.g. Goldsmith 
1998) and form an orthonormal set over the infinite Euclidean 
plane. The radial index can take the values p = 0,1,2, and 
the azimuthal index m can take the values m = 0, ±1, ±2, . . . ,±°°. 
In practice, it is customary to truncate the mode set by impos- 
ing maximum values for p and \m\. Indeed , it is usual to choose 
Hmax = Pmax, which leads to a total number of modes = 
(Pmax+ l)(2p m ax + 0- If /(r, 6) is real then clearly a p "' = (a»)*. 

In the context of quantum mechanics, it is more usual to label 
the modes in terms of the energy quantum number n = 2p+\m\ and 
m, rather than p and m. Indeed, Refregier (2003) uses this conven- 
tion to arrive at the polar shapelet modes 



CM) 



(2-5 0ffl )(q^)! 



1/2 



ImB 



(12) 

which also form an orthonormal set over the infinite Euclidean 
plane. 1 The energy quantum number can take the values n = 
0, 1,2, and the azimuthal (angular momentum) quantum 
number takes the values m = — n, —n + 2, . . . , n — 2,n, thereby giv- 
ing n + 1 values of m for each value of n. In practice, one truncates 
the mode set by imposing a maximum value for n, which leads to 
a total number of modes 5V2 = \ ("max + 1 ) ("max + 2) . Once again, 
if f(r, 0) is real, it is clear that b~ m = 

4.2 Properties of the mode sets 

It is clear from the above discussion (and that given in the Ap- 
pendix) that the two decompositions Jilt and 1121 use very dif- 
ferent combinations of Gauss-Laguerre modes. The question thus 
arises as to the relative merits of the two mode sets in decomposing 
a two-dimensional astronomical image. It is therefore of interest to 
investigate the properties of pixelised versions of these two mode 
sets. 

As shown by Refregier (2003), by analogy with the criteria 



1 We note that this definition of the polar shapelet functions differs some- 
what from that given in Refregier (2003) and Massey et al. (2003), which 
contain some typographical errors. 



£5), me appropriate values of the scale parameter (3 and the maxi- 
mum energy quantum number n max for the polar shapelet mode set 
<12> are given by 



P ~ (EWA 



vl/2 
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(13) 



for the decomposition of an image with features on scales ranging 
from mm to max . The galaxy images to be decomposed are of 
size 51x51 pixels, and contain structure down to scales of around 
3 pixels. Using the above criteria, we therefore choose rc m ax = 18. 
Thus, the total number of polar shapelet modes is 190. To allow a 
fair comparison between the polar shapelet and pm-mode sets, one 
should arrange for the two sets to contain (as close as possible) the 
same total number of modes. Fortunately, by choosing p m ax = 9, 
one can arrange for the pm-mode set also to contain 190 modes. 
The corresponding pm-mode and polar shapelet mode vectors are 
obtained by pixelising the continuous mode functions on the same 
grid as the galaxy images. Each mode set contains only 190 modes, 
which is clearly far less than the number of pixels in the image, but 
need not be less than the number of degrees of freedom in some 
class of images. For an arbitrary image of infinite resolution the 
number of degrees of freedom and pixels will be the same and both 
decompositions can only represent an approximation to the image. 

We may investigate the properties of the polar shapelet and 
pm-mode sets in the same manner as we analysed the one- 
dimensional shapelet modes. To this end, for each integer value 
of P in the range 0-20, we construct the mode matrix E for both 
the polar shapelets and pm-modes. The rank r of both matrices was 
found to equal the number of modes (190) for all values of p un- 
der consideration, hence showing that both modes sets are linearly- 
independent in all cases. To investigate whether the mode vectors 
in each set are orthonormal, we calculate the Gram matrix R=E'E 
for each mode set for each value of p. For an orthonormal mode set 
R should equal the identity matrix. In Fig. ll3l (top and middle), we 
plot the logarithm of the value of the largest and smallest diagonal 
element of R for the two mode sets. We see that, for both mode sets, 
the two curves coincide at log (/?,-,■) = only for P = 3 and 4, and 
so it is only for these values of the scale parameter that the mode 
vectors are of unit length. In the bottom panel of Fig. 1131 we plot 
the logarithm of the largest off-diagonal element of R for the polar 
shapelet (dashed line) and pm-mode set (solid line). We see that, 
for each mode set, the largest off-diagonal element is only below 
the machine precision for p = 3, and so it is only for this case that 
either mode set is orthogonal. 



4.3 Decomposition of HDF galaxy images 

For the galaxy images to be decomposed, the result 1 1 3i suggests 
one should set p = 14 for the polar shapelets; for comparison pur- 
poses we will also use this value of the scale parameter for the 
pm-mode set. 

To investigate the relative merits of the two mode sets, we 
first consider the decomposition of the single galaxy image shown 
in Fig. 1141 (left). In the middle panel, we plot the decomposition 
residual using the standard approach for the polar shapelets (dashed 
line) and the pm-modes (solid line) as a function of the total num- 
ber of modes M in the set. As expected, for both mode sets, the 
decomposition residual decreases as more modes are used. We also 
see, however, that, for the same total number of modes, the pm- 
modes yield a slightly lower decomposition residual than the polar 
shapelets. In the right-hand panel of Fig |14l plot the corresponding 
results obtained using the SVD appraoch. It is clear that, for any 
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Figure 13. Top: the logarithm of the value of the largest (solid line) and 
smallest (dashed line) diagonal elements of the Gram matrix R = E E as a 
function of the scale parameter p for the polar shapelet mode set. Middle: 
the same, but for the pm-mode set. Bottom: the logarithm of the largest 
off-diagonal element of if as a function of p for the polar shapelets (dashed 
line) and pm-mode set (solid line). 



given value of M, the SVD approach outperforms the standard ap- 
proach. We also see that, once again, the pm-modes yield a lower 
decomposition residual than the polar shapelets. 

The decomposition of five representative galaxy images into 
the pm-mode set are shown in Fig. l 151 using both the standard tech- 
nique (middle column) and the SVD method (right column). As 
anticipated, since the mode set is not orthogonal, the two decom- 
positions differ, markedly in some cases. For each galaxy image, 
however, we see that the SVD method produces a more faithful re- 
construction. In particular, we note that the SVD method is more 
successful in reproducing the finer detail of the image. By contrast, 
the decompositions obtained using the standard approach contain 
structure only on larger scales. A quantitative comparison of the 
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Table 1. The residuals e of the decompositions of the galaxy images shown 
in Fig |13l for the standard method and the SVD method. 



decomposition quality is given in Table Q which lists the values 
of the residual e of the decomposition. As expected, in each case 
the SVD method produces a more accurate decomposition than the 
standard method. 



5 DISCUSSION AND CONCLUSIONS 

We have discussed a general method for performing optimal im- 
age decomposition into a set of arbitrary digitised mode functions, 
which is an everyday problem in astronomy. Our approach may be 
applied straightforwardly even if the modes are linearly-dependent, 
non-orthogonal or incomplete, although in the last case the decom- 
position can only, in general, approximate the original image. By 
constructing the matrix £ containing the digitised modes, the prop- 
erties of the mode set can be easily obtained by performing a sin- 
gular value decomposition (SVD). Moreover, having obtained the 
SVD, the optimal values for the coefficients in the linear decompo- 
sition can be obtained immediately. If desired, the optimal coeffi- 
cients may also be obtained by constructing the set of dual modes, 
and performing scalar products of these duals with the input data 
vector. This should be contrasted with the standard method of tak- 
ing scalar products of the original mode vector with the data vector. 
This latter approach yields optimal results only in the case where 
the digitised mode vectors form an orthonormal set. In particular, 
we note that the digitisation process itself may lead to a mode set 
that does not inherit the properties of the continuous mode func- 
tions from which they are derived. Therefore, considerable care 
must be exercised when one performs decompositions using digi- 
tised versions of continous mode functions that form orthonormal 
sets. Adopting the standard approach without first investigating the 
properties of the digitised modes can lead to unnecessarily inaccu- 
rate image decompositions. 

We illustrate our general approach by first applying it to the 
decomposition of simple one-dimensional functions into Gauss- 
Hermite (shapelet) modes. We show that, although the continu- 
ous shapelet mode functions form an orthonormal set on the entire 
real line, the corresponding digitised mode vectors can be linearly- 
dependent and non-orthogonal, depending on the relative values of 
between the scale parameter p, the pixel size and the size of the 
image under analysis. We show that, for a wide range of these val- 
ues, the method developed here produces decompositions that are 
significantly superior to those of the standard approach. In some 
cases, the decomposition residual for our method is many orders of 
magnitude lower than that for the standard technique. 

We also illustrate our method by decomposing images 
of galaxies extracted from the Hubble Deep Fields into two- 
dimensional Gauss-Laguerre modes. We consider both the polar 
shapelet approach where the modes are indexed using the energy 
and angular momentum quantum numbers n and m, and the pm- 
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Figure 14. Left: a Hubble Deep Field galaxy. Middle: the decomposition residual using the standard approach for the polar shapelets (dashed line) and the 
pm-mode set (solid line) as a function of the total number of modes M in the set. Right: as middle panel, but using the SVD approach. 



modes used in optics, which are indexed by a radial number p and 
azimuthal number m. In both cases, although the set of continuous 
functions is orthonormal over the entire Euclidean plane, the corre- 
spondinig digitised modes can be severely non-orthogonal, depend- 
ing on the relative sizes of the mode functions, the image size and 
the pixel scale. We show that the SVD method proposed here again 
outperforms the standard approach, yielding more faithful decom- 
positions with lower residuals. We would therefore advise that in 
future applications of the shapelet decomposition method, the stan- 
dard approach should be replaced with that presented here. We also 
find that the pm-mode set provides consistently more accurate de- 
compositions than the polar shapelets. 

The decompositions presented in this paper used the SVD 
to provide an optimal decomposition for individual images. Fre- 
quently in astronomy, an analysis of the statistics of a collection of 
images is used to classify astronomical sources and properties. For 
such an application it may be desirable to use the duals method we 
have presented. In this case, the set of dual vectors can be calcu- 
lated using a SVD, without reference to any image. The expansion 
coefficients for each image can then be very efficiently calculated 
by performing scalar products with the duals. Additionally, at this 
stage, a priori knowledge can be used (if desired) to constrain the 
duals to preferentially extract certain features from any image. 

We conclude that the ability to work straightforwardly with 
linearly-dependent, non-orthogonal mode sets can be useful not just 
in accommodating digitisation effects. In many applications, it may 
be more efficient to work with modes that are known a priori to have 
such properties. For example, in one or two dimensions, it can be 
prove useful to decompose images simultaneously in shapelet mode 
sets centred at a number of points in the image. The 'multi-centre' 
shapelet bases will be discussed in a forthcoming paper. Finally 
we draw the readers attention to our own previous work investigat- 
ing generalised bases or arbitrary completeness and orthogonality 
within the field of THz optics (Berry et al. 2003; Withington et al. 
2002), including the description of second-order statistics associ- 
ated with partially-coherent fields. 
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Figure 15. Hubble Deep Field galaxies (left column) decomposed into polar 
Gauss-Laguerre pm-modes with p max = 9 and (3 = 14 using the standard 
method (middle column) and the SVD method (right column). 



12 R.H. Berry, M.P. Hobson & S. Withington 



REFERENCES 

Bertin E., Arnouts S., 1996, A&AS, 117, 393 

Golub G., van Loan C, 1996, Matrix Computations, John Hopkins Univer- 
sity Press, London 

Hobson M. P., Jones A. W., Lasenby A. N., 1999, MNRAS, 309, 125 

Massey R., Refregier A., Conselice C.J., Bacon D.J., 2003, MNRAS, in 
press I astro-ph/0301449 i 

Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P, 1994, Numeri- 
cal Recipies in Fortran 77, CUP, Cambridge 

Refregier A., 2003, MNRAS, 338, 35 

Refregier A., Bacon D., 2003, MNRAS, 338, 48 

Sanz J. L., Argiieso F., Cayon L., Martinez-Gonzalez E., Barreiro R. B., 

Toffolatti L., 1999, MNRAS, 309, 672 
Sanz J. L., Barreiro R. B., Cayon L., Martinez-Gonzalez E., Ruiz G. A., 

Diaz F. J., Argiieso E, Silk J., Toffolatti L., 1999, A&AS, 140, 99 
Tenorio L., Jaffe A. H., Hanany S., Lineweaver C. H., 1999, MNRAS, 310, 

823 

Williams R. et al., 1996, AJ, 112, 1335 
King I. R., 1972, AJ, 174, 123 

Murphy J. A., Withington S., 1996, Infrared Physics and Technology, 37, 
205 

Withington S., Murphy J.A., Yassin G, 2000, IEEE Trans. Antennas and 
Propagation, AP, 49 

Williams R. et al., 1998, A&AS, 193, 7501 

Berry R.H., Withington S., Hobson M.P., 2003, JOSA, in press. 

Withington S., Berry R.H., Hobson M.P., 2003, JOSA, in press. 

Berry R.H., Withington S., Hobson M.P., Yassin G, 2003, in Proc. of Four- 
teenth International Symposium on Space Terahertz Technology. 



APPENDIX A: PM-MODES AND POLAR SHAPELETS 

The Schrodinger equation for the two-dimensional isotropic har- 
monic oscillator may be written in plane polar coordinates as 

d 2 M/ 1 d\]f 
dr 2 r dr 

where we have set h = m = 1 and chosen our radial coordinate 

1 -2 



(13) (12) (11) (10) 



i a 2 \|/ 



- r V|/ = 2£\|/, 



such that the potential is given by V(r) = jr 1 . Seeking a separated 
solution of the form 6) = R(r)&(Q), one quickly finds the set 
of solutions 

\|^(r,6) =Aexp(-r 2 /2)(r 2 )l m l/ 2 4 m| (r 2 )exp(ime), 

where l)™\x) is an associated Laguerre polynomial and A is a nor- 
malisation constant. In addition to being energy eigenstates, these 
solutions are clearly also eigenstates of the angular momentum op- 
erator — i3/36 with eigenvalue m. The angular momentum quan- 
tum number may take the values m = 0, ±1, ±2, . . . , ±o°, whereas 
the 'radial' quantum number may take the values p = 0, 1,2, . . . , °°, 
Fig. lAll shows the mode function distribution in the pm-plane. 

Interestingly, the functions 1 12t are also the solutions of the 
paraxial wave equation in optics (see e.g. Goldsmith 1998), where 
they are called the Gaussian beam pm-modes or simply pm-modes. 
These modes are often used to describe the electric field distribution 
in paraxial quasi-optical systems. In this context, it is customary to 
truncate the mode set by imposing both a p max -value and a |m| max 
value. Moreover, it is usual to choose |m| max = p max , so the trun- 
cated set of modes form a rectangular in pm-space, as illustrated in 
Fig. lAll In this way, one obtains a total number of modes given by 

*C\ = (Pmax + l)(2pmax + !)• 

It is straightforward to show that the energy of the mode 
is given by (reinstating h for the moment) 

E = (2p+|m| + \)h = (n+l)h, 
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Figure Al. The distributions of mode functions in the pm-plane. The num- 
bers in brackets denote the value of the energy quantum number n for each 
mode. The dashed line indicates the usual geometry adopted in the pm- 
plane for truncating the mode set; in this C3.se pmax — I I max — 
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Figure A2. The distributions of mode functions in the ram-plane, The num- 
bers in brackets denote the value of the radial quantum number p for each 
mode and indeed may be a more intuitive index than n when investigat- 
ing images due to p being an intrinisic measure of radial extent and scale. 
The dashed line indicates the usual geometry adopted in the nm-plane for 
truncating the mode set; in this case n raax = 3. 



where n is the energy quantum number. Indeed, in the context of 
the two-dimensional harmonic oscillator, it is more usual to label 
the modes in terms of n and m, rather than p and m. In this way one 
arrives at the polar shapelet modes 

C(r,e)=Sexp(-r 2 /2)(r 2 )l m l/ 2 L^(r 2 )exp(ime), 

where B is a normalisation constant. The energy quantum number 
can take the values n = 0, 1 , 2, . . . , °°, whereas the azimuthal (angu- 
lar momentum) quantum number takes the values m = — n, —n + 
2, .. .,n — 2,n; hence the energy level with quantum number n is 
(« + l)-fold degenerate. Since we have performed just a simple 
relabelling, the polar shapelets 1121 and the pm-modes have the 
same functional forms, and are directly related by = <j>2^ + | m | • 
In Fig.lAH we show the distribution of modes in the nm-plane. In 
Refregier (2003) and Massey et al. (2003), to perform the decom- 
position of two-dimensional images, the mode set is truncated by 
imposing an rc max -value, as illustrated in Fig. lA2l In this way, one 
obtains a total number of modes given by 



5V2 = 5 ("max + 1 ) ("max + 2) . 
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